arXiv:astro-ph/9701126vl 17 Jan 1997 


Mon. Not. R. Astron. Soc. 000, jl] ^ (yr) Printed 1 February 2008 (MN IATj^X style file vl.4) 


Evolution of the Globular Cluster System in a Triaxial 
Galaxy 

R. Capuzzo-Dolcetta , 1 A. Tesseri , 1 

1 Istituto Astronomico, Universita La Sapienza, 

Via G. M. Lancisi 29, 1-00161, Roma, Italy 
dolcetta@astrmb.rm.astro.it, tesseri@astrmb.rm.astro.it 


Accepted yr m d. Received yr m d; in original form yr m d 


ABSTRACT 

Dynamical friction and tidal disruption are effective mechanisms of evolution of glob¬ 
ular cluster systems, especially in non-axysimmetric galaxies with a central compact 
nucleus. With a semi-analytical approach based on the knowledge of the dependence 
of the dynamical friction and tidal disruption effects on the relevant parameters, we 
are able to follow the time evolution of the globular cluster system in a model of a 
triaxial galaxy and give its observable properties to compare with observational data. 

An important result is that the flatter distribution of the globular cluster system 
relatively to that of the stellar bulge, as observed in many galaxies, can be explained 
by the evolution of the globular cluster system, starting from the same density profile. 
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1 INTRODUCTION 

Two observational facts are well established, by now: 

i ) the first is that many galaxies have globular clusters 
systems (GCSs) with density profiles less concentrated than 
their parent galaxy halo light, M87 and M49 being the pro¬ 
totypes (Lauer & Kormendy 1986, Harris 1986,1991). This 
has been recently confirmed by observations of a sample of 
14 elliptical galaxies, made with the WFPC2 of the Hubble 
Space Telescope (Forbes et al. 1996). These observations, 
thanks to the high resolution of the HST, were able to probe 
the inner kiloparsecs of those galaxies and show that most, 
if not all, of them have GCSs with surface density profiles 
that rise towards the centre less steeply than the underlying 
galaxy light; 

ii) the second is that many galaxies host compact mas¬ 
sive nuclei in their centres (Dressier & Richstone 1988, Kor- 
mendy 1988, Kormendy & Richstone 1995, Eckart & Genzel 
1996) with estimated masses in the range from 2-10 6 Mq for 
our Galaxy and M32 up to 3-10 9 Mq for M87. 

An explanation for the difference between halo star and 
globular cluster distributions has been proposed by Harris 
& Racine (1979), Harris (1986), and Racine (1991) as a dif¬ 
ference in the formation epoch of the two components. In 
this scenario, the GCS formed in an earlier phase of the 
protogalactic collapse, while the stars that constitute the 
halo condensed later, this resulting in a less peaked distri¬ 
bution for the clusters. This picture requires an exact tim¬ 
ing in the sequence of the evolution, in order to permit the 
clusters to be less metal rich than the halo stars, while pro¬ 


ducing the required differences in central concentration. In 
disk galaxies,however, ’the epoch of cluster formation would 
be early enough to force chenical enrichment but not early 
enough to take on a distinct spatial structure’ (Harris 1986, 
Sect. VII, p. 840). Moreover, this scenario does not explain 
why the tails of the two density distributions are almost the 
same. This last observational evidence suggests an alterna¬ 
tive explanation: the cluster system and the halo formed 
at the same time with a similar spatial distribution, and 
the present differences are a consequence of the dynamical 
evolution of the GCS. Dynamical evolution correlates also 
with the possible presence of masssive central nuclei. Ac¬ 
tually, the larger core radius of the GCS would imply that 
the globular cluster population has been significantly de¬ 
pauperated in the inner regions of the system. This is the 
case when a massive object (like a black hole) resides in the 
centre of a galaxy and disrupts, by means of tidal forces, 
globular clusters which pass sufficiently close to it. There is 
no direct evidence that the aforementioned massive objects 
are black holes (they could be, as proposed by Kormendy 
& Richstone (1995), massive clusters of low-mass stars, stel¬ 
lar remnants etc.) except in the case of NGC 4258, where 
the discovery of a perfect keplerian rotation curve in the 
inner regions (Miyoshi et al. 1995) rules out, on dynamical 
evidences (Mayoz 1995), alternatives to black holes. In the 
outer regions, on the contrary, the influence of the central 
massive object will be negligible. So we expect that there the 
cluster distribution has remained more or less unchanged. 

Another major dynamical effect, the dynamical friction 
of field stars on globular clusters, enhances the efficiency of 
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the depleting mechanism. It acts reducing the cluster or¬ 
bital energy and bringing the clusters towards the centre, 
thus increasing the number of globular clusters in the inner 
regions. These clusters could feed the massive object. This 
scenario has been proposed first by Tremaine, Ostriker and 
Spitzer (1975), in a study on M31, where they showed how 
a massive object of mass in the range 10 7 -10 8 Mq could 
directly form from globular clusters braked to the centre of 
the galaxy and there merged. 

Both of these effects (dynamical friction and tidal dis¬ 
ruption) are significantly emphasized if the galaxy is triax- 
ial in shape, a possibility supported by several observations 
(see, for example, Bertola, Vietri & Zeilinger (1991) which 
show evidence of triaxial distributions in 32 galaxies). The 
orbits which constitute the bulk of such a potential are the 
’box’ orbits, which are dense around the centre (see, e. g., 
Binney & Tremaine 1987, hereafter BT) where the massive 
object lies (so that even globular clusters of large apocen- 
tric distance will possibly be disrupted) and the field star 
density is higher (so that dynamical friction efficiency is in¬ 
creased). In fact, Pesce, Capuzzo Dolcetta and Vietri (1992) 
have demonstrated that dynamical friction decay times on 
box orbits are significantly reduced. 

The shape of the velocity ellipsoid of halo stars and 
globular clusters in our galaxy supports this picture. In fact, 
while the velocity dispersion of the halo stars is larger in the 
radial direction, as expected from numerical simulation of 
the radial collapse of the protogalaxy, the globular clusters’ 
velocity ellipsoid is almost spherical. Under the hypothesis of 
a coeval formation, this can be explained by a selective pro¬ 
cess which destroyed the clusters on low-angular momentum 
orbits (i.e. box orbits in a triaxial galaxy). The role of triax- 
iality in the tidal disruption mechanism has been quantita¬ 
tively discussed by Ostriker, Binney and Saha (1989) (here¬ 
after OBS). In their paper they proposed, for the first time, 
that the formation of a massive nucleus from decayed globu¬ 
lar clusters could be a self-limiting process, due to the inverse 
proportionality of the tidal disruption timescale, Tud, to 
the nucleus mass. Capuzzo Dolcetta (1993) has investigated 
thoroughly the evolution of a GCS in the Schwarzschild’s 
(1979) triaxial non-rotating model, under the combined ef¬ 
fects of dynamical friction and tidal disruption, studying the 
growth of the nucleus and the evolution of globular clus¬ 
ters’ mass function. Indeed, he found that the cooperation 
of these two effects may lead to the formation of a compact 
nucleus, in form of globular clusters decayed to the centre of 
the galaxy. The growth of the nucleus eventually halts when 
its mass is large enough to shatter all incoming clusters. Of 
course the value of the mass reached by the growing nucleus 
depends, in this scheme, on the initial GCS spatial, mass 
and velocity distributions. 

The increased efficiency of dynamical friction and tidal 
disruption, in a triaxial galaxy, is so increased that their 
effects are not limited to the very inner regions of the par¬ 
ent galaxy. Moreover, there is no need to assume, for the 
GCS, a box-biased phase-space density such that the globu¬ 
lar clusters are all on box orbits. In fact, Capuzzo Dolcetta 
(1993) showed that even in the case of an isotropic distri¬ 
bution function (hereafter DF), the evolution of the GCS is 
very similar to that of a box-biased DF. The reason for this 
is that the two DFs do not differ much in the region of the 
phase-space occupied by the majority of the clusters. 


In Section 2 we describe our model of a globular cluster 
system evolving due to dynamical friction and tidal disrup¬ 
tion effects and we give a formula which permits to calculate 
the density profile of the GCS and its observable properties; 
in Section 3 we present and discuss the results. 


2 THE MODEL 


In this paper we study a population of clusters exclusively 
on box orbits. This may be the case if they formed in the 
early galactic stages, during the radial collapse of the proto- 
galactic nebula (van Albada 1982, Binney 1988). We develop 
a semi-analytical model which allows to follow the evolution 
of the spatial density of the GCS in a triaxial galaxy under 
the influence of the two main evolutionary effects: dynamical 
friction and tidal disruption. Once the spatial density profile 
of the GCS is obtained, we can deduce its surface density 
profile, core radius and other useful quantities which may 
be compared to observations. The galactic potential adopted 
here, like in Pesce et al. (1992) and Capuzzo-Dolcetta (1993), 
is that of the Schwarzschild’s (1979) model. We define (fol¬ 
lowing Capuzzo Dolcetta 1993) Tdf(E,m) as the time re¬ 
quired to a cluster on a box orbit to lose all its energy E 
(and stop at the centre of the galaxy) by means of the fric¬ 
tional drag exerted by the stellar population. A good fit to 
Tdf(E,m) is: 

T df (E, rn ) = - Tr - w yr, (1) 

where m is the mass of the cluster in units of 10 6 Mq and E 
(0< E < 1) is the orbital energy per unit mass in units of $o, 
the central value of the Schwarzschild potential. Equation 
(1) is not well behaved at the origin since, as the energy 
approaches zero, Tdf(E,m) reaches a Unit limit, which is 
inconsistent with its definition. In practical calculations we 
used a slightly different form of Eq. (1), which exhibits the 
behavior required. 

To include the effect of tidal disruption we will need 
an estimate of the timescale, Tud, associated to this effect 
which, according to OBS and Capuzzo Dolcetta (1993), may 
be expressed in the form 


Ttid, — 


Gm v n A v 


fj,7T\/5 V R h vl r c R h 


Tr 


( 2 ) 


where: u = GA h ■ M n is the mass of the nucleus; r c is the 

^ r c v£ 

core radius of the galaxy; v c is the circular speed at large r; 
v n is the speed of the cluster of mass m at the point where 
the gravitational attraction of the ellipsoid equals that of 
the nucleus; A w is the area of the waist of the box orbit; Rh 
is the half-mass radius of the cluster; T r is the half-period 
of oscillation parallel to the potential long axis. Note that 
Eq. (2) depends on the structural parameters of the clusters 
only through the quantity /-S- which is proportional to 

V H h 

the square root of the mean density within the half-mass 
radius, \fph- The dependence of Tud on the orbital energy 
E is through A w , v n and T r . 

Since we are interested in the density profiles of the 
GCS, it is useful to treat it as a ‘fluid’ system. This allows 
us to adopt the continuity equation as the one which rules 
the evolution of the system: 
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dp 8pV r pv r 
dt + dr + r 


(3) 


where p(r, t) is the number density of the GCS, v r {r, t ) is its 
radial velocity field and S > 0 is a ‘sink’ term. We have writ¬ 
ten the continuity equation in spherical coordinates, drop¬ 
ping the terms containing the derivatives with respect to 
angles, since the observed GCSs usually show spherical sym¬ 
metry. Moreover, we have adopted vg =0 and v^=0 since the 
GCSs seem not to rotate. If no evolution of the GCS occurs 
(as is the case when the system is collisionless and no mas¬ 
sive central ‘absorber’ is present), we expect Ur-(r,t)=0 and 
S=0 for any r and t. In our picture, instead, the frictional 
drag of the stellar halo on the GCS, reduces the energy of 
the clusters, acting as a radial velocity field pointing inward. 
At the same time, the massive galactic nucleus erodes the 
GCS, standing for a non-zero sink term S = p/rud, where, 
hereafter, Tud, is the average of formula (2) over the DF 
of the system, as described in Appendix A. Once we give 
the initial GCS density distribution p(r, 0), and we assume 
v r (r,t)=v r (r) (that is, the radial velocity field generated by 
the unevolving stellar halo does not change during the evo¬ 
lution of the system), the following solution of Eq. (3) is 
obtained (see Appendix B): 


p(r,t) = 


X 


where 


v r (H ( r,t))H^ ( r,t ) 
v r (r)r a 


p(H(r,t), 0) x 



_ dx _ 

T t id( H ( r ,x),M n (t-x),p h ) 


(4) 


H(r,t) = H (t + r{r)) (5) 

T(r) '"/w) (6) 

and H(x) is the inverse function of r(r). The function M n (t), 
which gives the nucleus mass at any time t, is not known. 
In our calculations, we will adopt a fixed nucleus mass since 
a reliable evaluation of the rate of accretion of the nucleus 
requires further investigation (see Capuzzo-Dolcetta 1996). 
Anyway, when a realistic model will be available to give 
Mn(t), it will be straightforwardly included in our model. 

By the definition of r(r), the difference t(R) - r( 0) is 
the time required to an element of mass of the ‘fluid’ to reach 
the centre and stop there, starting from a point at distance 
R from the origin. It is so straightforward to relate this r(r) 
to the function Tdf( E,m). This is easily done by averaging 
Tdf ( E , m) over the DF of the system, thus obtaining a func¬ 
tion Tdf( r,m) (see Appendix B). Now, the velocity field at 
any point r can be estimated by 

T 

v r {r,m) = -. (7) 

T df 

Thus, substituting Eq. (7) in Eq. (6), we find the following 
expression for the function r(r) which appears in the solution 
(4): 


r(r, m) = J -^-dr. (8) 

In this way we have related the dynamical friction effect to 
the presence of a velocity field pointing radially inward. As 
expected, p(r,t) as given by Eq. (4), depends on m via Eq. 
(8). In the following, we assume the initial distribution of 
the GCS in the form 


© yr RAS, MNRAS 000, | § 


p(r,Q;m) =tp 0 (m) • p(r,Q) (9) 

where m is the individual cluster mass and i/>o(m) is the 
GCS’s initial mass function (hereafter IMF). The relation 
(9) says that the IMF of the GCS is indipendent of the 
position in the system. The solution of eq. (3), in this case, 
is given by a superposition of distributions: 

j-m 2 

p(r,t) = / p(r,t;m)ip(m)dm (10) 

J m\ 

where p(r,t;m) is the solution (4) for a single mass GCS, 
taking into account the dependence of Tdf on the clusters’ 
mass, m. To conclude, the dependence on the overall model 
is given by the functions Tdf and Tud, and by the GCS’s DF 
which fixes its initial distribution. It is clear, hence, that 
equations (4) and (10) may represent the evolution of a GCS 
in different situations, provided we know Tdf, Tud and f(E). 


3 THE RESULTS 

In the previous Section we have solved the evolution equa¬ 
tion for the density distribution of a GCS made up of clusters 
of different masses when dynamical friction and tidal disrup¬ 
tion are taken into account. Before showing the results, we 
specify the parameters which our model depends on. 

For ph, we seek a relation with the cluster mass. From 
data of clusters in our Galaxy (Webbink 1985) we found 
a loose relation among the mass and the half-mass den¬ 
sity of the clusters, obtained by adopting a fixed ratio 
(M/L v )q= 1.5 and fitting the clusters density profiles with 
Plummer models: 

p h = 1.38 • 10 5 m 1 ' 56 M e pc~ 3 (11) 

This choice for the cluster half-mass density is quite dif¬ 
ferent from that used by OBS in their calculations. Their 
GCS is composed of identical clusters with half-mass den¬ 
sity ph = 10 3 Mgpc -3 , implying that we are dealing with 
clusters significantly denser than theirs. Since the tidal dis¬ 
ruption timescale depends on yfph, we expect a lesser deple¬ 
tion of the GCS. 

As initial single-mass number density for the GCS we 
assume: 

p(r,°) = --£2--y (12) 

1 + ( —) 2 2 

which corresponds to the monopole component of the 
Schwarzschild potential, po being the central density and 
r c o the initial core radius (assumed equal to the bulge star 
core. Note that the projection, E(r), of (12) is the so-called 
modified Hubble law, a function well fitting the elliptical 
galaxies’ surface brightness: 

E(r ’ 0) "77TZy (13) 

' V ICO / 

where Eo = 2por c o- Since the density distribution (12) has 
an infinite mass we cut it at a radius r max 3>r c o. Varying the 
value of l'max results in a different behaviour of the averaged 
timescales Tdf and Tud■ Anyway, we find that, letting l'max to 
vary in a reasonable range, the differences are more impor¬ 
tant in the outer regions ( r > 20 r c o), where, anyway, the 
timescales are always sensibly longer than a Hubble time. 
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Figure 1. Evolution of the GCS space and surface density (upper 
and lower panels,respectively) due to dynamical friction only. Left 
panels: flat (s=0) IMF; right panels: steep (s=2) IMF. The thick 
curve is the initial distribution, the other curves correspond to 1, 
5, 10 and 15 Gyr (bottom to top) 
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Figure 2. Same as Fig.l, taking into account tidal disruption of 
clusters by a nucleus of fixed mass M rl =10 7 Mg. The thick curve 
is the initial distribution, the curve corresponding to 1 Gyr is the 
closest in shape to the initial distribution and the other curves 
correspond to 5, 10 and 15 Gyr (top to bottom) 


In the regions of interest (r< 10 r c o), the changes in r max 
reflect in negligible changes in our timescales. 

As DF for the GCS, we use a King model with <j 2 =8 ( I>o 
(see eq. 4-130 and 4-131 in BT), whose corresponding den¬ 
sity profile fits well eq. (12), as required. 

For the GCS IMF we shall assume truncated power- 
laws: 

! 0 0 < m < mi 

k m~ s mi < m < m 2 (14) 

0 m > m 2 

with mi = 10 4 Mq and m 2 = 3- 10 6 Mq, so that the product 
ipo{m)p(r, 0) gives the intial number per unit volume of clus¬ 
ters with mass m. We consider the two cases of a flat (s=0) 
and a steeply decreasing (s=2) IMF. The normalization fac¬ 
tor kpo is chosen in such a way to give a total number of 
clusters lVt o t=1000. 

3.1 Dynamical friction 

When a pre-existing nucleus is absent - or its mass is so 
small that the tidal disruption effect is negligible - we may 
apply formula (B.ll) for the evolution of the number density 
distribution of a GCS undergoing dynamical friction only. 
The results are shown in Fig. 1. We observe an inner re¬ 
gion (within r« r c o)where the density is strongly enhanced, 
surrounded by a depleted strip (r c o < r < 10 r c o). This ef¬ 
fect is best seen in the s=0 case (Fig. la,c), since, with such 
an IMF, the GCS is composed by a large fraction of mas¬ 
sive clusters, which decay faster. On the contrary, with a 
steep (s=2) IMF (Fig. lb,d), the GCS is composed mainly 
of light clusters, which respond slowly to the frictional drag, 


thus the depletion in the outer region is almost impercepti¬ 
ble. The surface distributions, in both cases, become more 
concentrated as time goes on, so that the core radii become 
smaller. This is a natural consequence of the fact that globu¬ 
lar clusters lose their energy by dynamical friction and move 
on less extended orbits. 

To study the time evolution of the core radius of the 
GCS, we need to remember that an actual observation can¬ 
not sample the GCS distribution all the way to the centre 
of the galaxy. As a consequence, the value of the observed 
core radius depends on a somewhat arbitrary extension of 
the GCS’ surface density towards the centre of the galaxy. 
Were the inner limit of the observations equal to, say, r c o, in 
the s=2 case we would not observe any differences between 
the halo and GCS distributions (Fig. Id); in the s=0 case we 
would observe a flatter slope (Fig. lc), that would lead us 
to conclude erroneously that the GCS is less concentrated 
than the halo. Could we sample clusters to a smaller inner 
radius, say 0.1 r c o, we would observe that the GCS istribu- 
tions, in both cases, are more concentrated than that ofthe 
halo. It is evident, then, the importance of going with the 
observations as close to the centre of the galaxy as possible. 
In the next Subsection we give a more detailed evaluation 
of the core radius of the evolved GCS. 

It is interesting to calculate the number of clusters con¬ 
tained in a given radius R at time t, N c i(R, t). In the case of 
dynamical friction alone, the integral over r is easily done, 
and we obtain: 

pn 2 

N c i(R, t) = / N°i(H (t + Tdf(R,m)))4’o(m)dm (15) 
J mi 

where N°i (r) is the number of clusters in the sphere of radius 
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Table 1. Number and total mass of clusters inside one core ra¬ 
dius at time 15 Gyr, compared to their initial values (with the 
subscript 0). 


IMF 

N c io 

N c z 

M c zo (Mg) 

M c z (Mg) 

s=0 

37 

410 

5.610 7 

7.2-10® 

s=2 

37 

66 

2.1-10 6 

1.2-10 7 
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Figure 3. Same as Fig.2, but M n =10 8 Mg. 

r at t = 0. A similar equation holds for the mass contained 
in the sphere of radius R at time t, M c i(R,t). In this way, 
we can evaluate the number and mass of clusters within the 
core at present time (f=15 Gyr), N c i and M c i, and compare 
them with their initial values, N c io and M c zo (see Table 1). 
This table gives a first insight on the possible formation of a 
massive nucleus by globular clusters accretion, since the M c ; 
values constitute the upper limits to the mass contribution 
by the cluster system. It is easily seen that N c i(R , f), as well 
as M c i(R,t), scale linearly with Ntot- So, bigger values for 
the mass inside the core radius require larger Ntot, which 
is not a free parameter, being constrained by observations, 
through the presently observed number of clusters. 

3.2 Dynamical friction and tidal disruption 

Now, we show the results for the volume density and surface 
distribution of a GCS subjected to both dynamical friction 
and tidal disruption. As it is easilyseen, in Figs. 2, 3 and 4 
the density profiles strongly depend on M„ 

The role of the nucleus is overwhelming when the GCS 
is composed mainly of low-density clusters (s=2 case, see 
Eq. (11)). In this case, the profiles differ significantly from 
the case of the absence of a nucleus, even if the nucleus mass 
is moderate (M n = 10' Mg, compare Fig. lb,d to 2b,d). On 
the contrary, the more massive clusters of the flat IMF, be¬ 
ing denser, resist effectively to the tidal disruption when 
M n < 10' Mg (compare Fig. la,c to 2a,c). Of course the de¬ 
le) yr RAS, MNRAS 000, | | 
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Figure 4. Same as Fig.2, but M n = 10 9 Mg. 

pauperating effect of tidal disruption is enhanced by heavier 
nucleus masses (see Fig. 3 and 4). 

Another interesting feature is the prominent central 
density cusp, more evident for light nuclei and flat IMF. 
This is due essentially to heavy clusters which rapidly move 
towards the centre of the galaxy due to large dynamical fric¬ 
tion suffered by their high mass, and contemporarily survive 
tidal disruption because of their high density. At increasing 
nucleus masses, this effect is less evident and becomes almost 
imperceptible for heavy nuclei, because of the increasing ef¬ 
ficiency of tidal disruption. Note how the differences among 
the density profiles of the innermost and external regions 
are significantly reduced by projection (compare upper and 
lower panels in Fig. 2,3 and 4). 

3.2.1 The GCS core radius evolution 

In the s=2 case for the IMF, the surface distributions de¬ 
part from a modified Hubble law just in the inner regions, 
thus leading to a reliable r c determination. With a flat IMF, 
instead, no region of constant density (a ‘core’) is kept up to 
present time. Anyway, observationally, the steepening of the 
profile (which carries the signature of dynamical friction) 
could be appreciated just inside a region whose radius is of 
the order of r c o- In this case, the evaluation of r c depends 
on the extension of the observed profile towards the centre, 
which usually results in its overestimate. To give a quantita¬ 
tive evaluation of this overestimate, we have calculated the 
’observable’ core radii of the GCS after a dynamical evo¬ 
lution of 15 Gyr, by fitting the surface density profiles of 
Fig. 5 with modified Hubble laws. By ’observable’ we mean 
that the fits were done excluding too inner regions, within 
some radius r m i n - The results are given in Tables 2 and 3, 
for different nucleus masses. In the s=2 case the value of r c 
does not change much on varying r m i n , confirming that the 
surface distribution of the GCS is well approximated by a 
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Figure 5. Plot of the GCS surface density profile after evolution 
has occurred for a flat (s=0) IMF (upper panel) and for a steep 
(s=2) IMF (lower panel). The thick curve is the initial profile, 
the other curves correspond to nucleus masses of 10 7 , 10® and 
10 9 Mg (top to bottom). 


Table 2. ‘Observed’ core radius for the case of a flat IMF (s=0). 


r min 

M n =10 7 

M n =10 8 

M„=10 9 

fco/2 

2.7 

8.7 

24.0 

r c o 

3.0 

9.3 

25.0 

5 r c0 

3.0 

9.3 

25.0 


Values obtained by fitting the surface density profiles of the 
evolved GCS with a modified Hubble law, excluding the region 
inside r m i n , for different nucleus masses M n . 

modified Hubble law. On the contrary, in the s=0 case, the 
observed core radius depends strongly on the value of r m i n , 
indicating that the surface distribution does not display a 
well defined central core. Then, a core is reliably defined 
when clusters have been destroyed by tidal disruption. 

In Fig. 6, we plot, for the steep IMF, the time evolution 
of the ‘observed’ core radius for different nucleus masses, 
obtained as described above. 

3.2.2 A parameter reliably comparable with observations 

A better parameter to synthesize the properties of the 
evolved GCSs is given by x = Ni/No, defined as the ra¬ 
tio between the number of ‘lost’ clusters outside the min- 

Table 3. Same as Table 2, but for the case of a steep IMF (s=2). 
r min M n = 10 7 M n =10 8 M„=10 9 


r c o/2 

0.8 

0.8 

1.6 

r c o 

1.4 

2.0 

4.2 

5 r c o 

4.6 

5.2 

7.5 



Figure 6. Time evolution of the ‘observed’ core radius for a GCS 
with a steep (s=2) IMF and M n = 10 7 Mq (dot-dashed curve), 
M n = 10 8 Mq (dotted curve), M n = 10 9 Mq (solid curve) 


Table 4. Fraction of ‘lost’ clusters after 15 Gyr. 


M n (Mq) 

s=0 

s=2 

10 7 

0.40 

0.20 

10 8 

0.42 

0.46 

10 9 

0.51 

0.78 


imum radius reachable by the observations (r min) to the 
initial number of clusters in that region. The number Ni is 
given by the difference between the initial and the present 
number of clusters, outside r m i n . The ‘lost’ clusters include 
either clusters which have been actually destroyed by tidal 
interaction with the massive nucleus either clusters which 
have lost energy by dynamical friction and have moved on 
orbits all within r m in. 

The quantity a: is a good parameter to quantify the roles 
of the evolutionary mechanism under discussion, because it 
could be easily inferred from observations, as the difference 
observed between the (normalized) radial profiles of halo 
stars and GCS. This requires the reasonable hypothesis that 
the initial distributions of clusters and bulge stars were the 
same and that the present bulge distribution is equal to the 
initial (see McLaughlin 1995, Capuzzo-Dolcetta & Vignola 
1996). In Table 4 we give the fraction of ‘lost’ clusters for 
a GCS subjected to dynamical friction and tidal disruption 
with different nucleus masses and for the two IMFs, with 
r m in = r c o- For the flat IMF, x varies of just 25 per cent at 
varying the nucleus mass from 10 7 Mq to 10 9 Mq, indicat¬ 
ing that dynamical friction is important. For the less dense 
clusters of the steep IMF x increases by a factor four over 
the same M n range because, in this case, tidal disruption is 
dominant. 

To understand better the contributions of the two de¬ 
pleting mechanisms and their dependence upon the relevant 
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parameters (nucleus mass and clusters’ half mass density) 
we have followed the evolution of the GCS density profiles 
in two limiting cases: i) pure dynamical friction and ii ) pure 
tidal disruption. The computations were performed on var¬ 
ious single-mass GCSs, over a range of individual cluster 
masses m, choosing = rco¬ 
in case i), the results at time 15 Gyr are well fitted by 
a broken power-law, 

f 0.46-m 0 ' 82 0.01 < m <0.6 

Xdf = \ 0.35-m 0 36 0.6 < m <3 (16) 

where m is in 10 6 A/q. For a GCS with a distribution of 
masses it is quite natural to calculate the quantity Xdf by 
averaging it on the IMF. This leads to an error of about 3 
per cent respect to the exact result obtained by integrating 
the volume density given by Eq. (10), indicating as reliable 
the simple average of the expression (16) on different IMFs. 
For our models, we obtain: for a flat IMF, Xdf= 0.38; for 
a steep IMF, Xdf= 0.04. Clearly the influence of dynamical 
friction is almost negligible if the GCS’s initial mass function 
is biased towards light clusters. 

In case ii), we expect that Xtu depends on the nucleus 
mass and cluster half-mass density through the combination 
Mn/sfph- In fact, we find that a good fit to Xtid is 

_ / 6.79 • 10 _5 y°’ 79 0 < y < 7.47 ■ 10 5 . . 

Xud ~f 0.17 ■ y 0 ' 21 y> 7.47-10 s ' /! 

where y = M n ■ ppf^ 2 , ph is in units of Mq- pc -3 and M n in 
solar masses. Like before, to obtain the quantity Xtid for a 
GCS with a distribution of masses we average the expression 
(17) on the IMF. 

When both the effects are at work, a reasonable expres¬ 
sion for x would be, clearly: 

x = Xdf + (1 - Xdf) xud- (18) 


The values obtained in this way are in good agreement with 
the s=2 case of Table 4, since Xdf <C xud, while the s=0 case 
is poorly represented because the two effects compete. 


3.3 The evolution of the mass function 


Since both dynamical friction and tidal disruption timescales 
depend differently on the cluster mass, the IMF given by 
Eq. (13) evolves into a mass function which is different from 
point to point in the system. 

In Fig. 7 we show the present time mass function of 
the GCS for different nucleus masses (excluding the region 
within r c o to make it comparable to radially limited observed 
samples). 

Note that, as long as the GCS mass contribution to the 
nucleus is negligible, the knowledge of the evolution of a 
particular mass function allows to obtain the evolution of 
a different IMF. Indeed, under this hypothesis, the individ¬ 
ual mass components of the GCS evolve independently of 
each other. Hence, the following relation among two mass 
functions ip and <p 


Ipe _ <Pe 
ipo <po 


(19) 


holds, where the subscript ‘0’ stands for the IMF and the 
subscript ‘e’ stands for the evolved mass function. Thus, 
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Figure 7. Evolved (t=15 Gyr) mass functions for different nu¬ 
cleus masses. Upper panel: steep (s=2) IMF; lower panel: flat 
(s=0) IMF. The thick solid curve is the IMF, the other curves 
correspond to M n = 10 7 , 10® and 10 9 Mq (top to bottom). Clus¬ 
ter masses are in 10 6 Mq; ip in units such that A7 fJ / = 1000. 


it is straightforward to obtain the initial mass function (po 
from the evolved observed mass function <p B via Eq. (19). 
Of course, some assumptions on the time evolution of the 
nucleus mass are needed. In a more detailed scenario, the 
GCS will contribute to M n in a way dependent on the IMF, 
which means that Eq. (19) no longer applies. 


4 CONCLUSIONS 

In this paper we have developed a model which allows to 
follow the evolution of the density distribution of a globular 
cluster system (GCS) in a triaxial galaxy under the influence 
of two effects: dynamical friction by field stars and tidal 
disruption by a massive central object. Both these effects are 
amplified by the triaxiality of the gravitational potential of 
the parent galaxy. The exact knowledge of the present-time 
density profiles permits to calculate the value of some of the 
relevant observables as a function of the assumptions made 
on the initial distribution and mass function of the GCSs. 

Actually, we find that the minor concentration of the 
GCSs, relatively to the distribution of halo stars, as ob¬ 
served in many galaxies, is an effect which depends strongly 
on the quality of the observations. In particular, a compar¬ 
ison of the core radius of the GCS with that (r c o) of its 
parent galaxy star distribution might be misleading, since 
the estimate for the core radius of the GCS depends on the 
minimum radius at which clusters are sampled. For example, 
we have shown that, in the case of a GCS whose evolution is 
prevalently ruled by dynamical friction, as it is the case for 
a GCS made up mostly of heavy clusters, the final shape of 
the distribution displays a strong concentration of massive 
clusters in the very inner regions of the galaxy and a con- 
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sequent lack of clusters in the outer regions. In this case, if 
the observations do not reach the inner regions (within r c o), 
we effectively measure an erroneusly large core radius. Thus, 
an observed larger core radius is not a firm signature of the 
presence of a massive object at the centre of the galaxy. 
However, we found noticeable differences among the case of 
a steep IMF and that of a flat IMF. These differences, due 
to the influence of the half-mass density of the clusters as a 
function of the cluster mass, are such that: 

i ) for a steep IMF, made up mostly of light clusters, 
the prevailing effect is tidal disruption and there is a strong 
dependence of the evolved core radius on the mass of the 
central object (the value r c = 10r c o is obtained with a nu¬ 
cleus of 2-10 8 Afg); 

ii) for a flat IMF, made up mostly of heavy (and dense) 
clusters, the prevailing cause of depletion is dynamical fric¬ 
tion and, consequently, the dependence on the mass of the 
central object is weak (r c = 10r c o for 2 • 10 9 Mq). 

We have found that a parameter which is more reliable 
than the core radius to describe the GCS is the number 
of ‘lost’ clusters outside some minimum radius (we choose 
r-min = r c o). This is the difference between the initial and 
the observed number of clusters, under the assumption that 
the GCS was initially distributed as the parent halo light. 
Our calculations show that: for a flat IMF the percentage 
of ‘lost’ clusters ranges from 40 per cent (with no massive 
central object) to 50 per cent (for a 10 9 Mq nucleus) of the 
initial total number; for a steep IMF it ranges from 3 per 
cent (no nucleus) to 80 per cent (10 9 Mq nucleus). 

We give two formulas which fit the fraction of ‘lost’ clus¬ 
ters in the case of dynamical friction only and in the case 
of tidal disruption only, as a function of the cluster mass. 
When both the effects are working, the number of ‘lost’ clus¬ 
ters may not be obtained by simply summing the fractions, 
confirming that an interaction among the two effects exists. 

To conclude, our calculations show that the present dif¬ 
ferences observed between the GCS and halo stars surface 
distributions can be explained, by dynamical evolution of the 
GCS , under the influence of dynamical friction and tidal dis¬ 
ruption, even if the initial concentrations (core radii) were 
the same. Thus, to have a definite answer to the question 
’’are the observed differences between the star-bulge and the 
GCS density profiles just reflecting different initial condi¬ 
tions or are a consequence of evolution?”, it would be crucial 
to compare their kinematical properties. For example, the 
knowledge of the run with radius of the GCS (projected) 
velocity dispersion may help to understand if, like in our 
galaxy, there has been a selective depauperation of clusters 
on highly radial orbits. 
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APPENDIX A: 

The relation 

r a = u~ 1 (E) (Al) 

holds between the orbital energy E and apocentric distance 
r a , where u(r) is the spherical symmetric component of the 
Schwarzschild’s potential (see Eq. (5) in Pesce et al. 1992). 
Eq. (Al) suffices to give Tud as a function of E. At any point 
r, we average the tidal disruption timescale on the distribu¬ 
tion function, f(E), of the clusters system, thus obtaining 
an averaged tidal disruption timescale: 

'Ttid{r', M n , Ph) — (xtid^E , M n , phf) DF — 

J T t id (E,M n ,p h )f(E) v 2 dv ( A2 ) 

j f(E) v 2 dv 

where the integration is done on the region of phase-space 
corresponding to bound orbits. In the same way, we may 
obtain an averaged dynamical friction timescale 

Tdf(r,m ) = (T d f(E)) DF = 

J T df (E)f(E) v 2 dv ( A3 ) 

J f (E) v 2 dv 
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APPENDIX B: 


First of all, we seek a solution, say p of the continuity equa¬ 
tion in its homogeneous form, that is: 


dp dpv r pVr _ 
dt + dr + r 


(Bl) 


An exact analytical solution is easily found when v r (r,t) = 
u(r)w(t), in the following way. 

We define the quantity N(r, t) as 


nr 

N(r,t)= / Ln / p{r' ,t)r' 2 dr 1 

Jo 


(B2) 
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Now, integrating Eq. (Bl) over the same spherical volume 
as in Eq. (B2), we obtain an equation for N(r,t ): 


dN dN 
~dt + 


= 0 . 


If we define the functions 


(B3) 



Eq. (B3) becomes 

dN dN_ 
dW dr ~ ' 
whose solution is: 


(B4) 


(B5) 


(B6) 


N(r,t) = F(W(t) + r(r)) (B7) 

where F(x) is a function constrained by the initial conditions 
on p(r,t). The solution /5(r, t) is obtained by inverting Eq. 
(B2): 

« r ' t >= 4 ^'^ = ^ /(H ' < ‘ )+T(r)) <B8) 

Assuming an initial density distribution p(r, 0), we find the 
relation 


f(x)=p(H(x),0)u(H(x))H(x) 2 (B9) 


where H(x)= r 1 (x). Finally, the solution of the homoge¬ 
neous continuity equation is 


p(r,t) = 


u(H(W(t) + r(r)))H 2 (W(t) + r(r)) 
u(r)r 2 


x p(H(W(t) + T(r)),0) 


(BIO) 


For the purposes of this paper, we set w(t)=l (see Sect. 2), 
so that W(t)=t and the solution (BIO) simplifies to: 


P{r,t) 


u.(H(t+T(r)))H 2 (t+T(r)) w 
- 7 — \ 7 - '*■ 

uyr)r 


X p(H(f + r(r)),0) 


(BH) 


When a ‘sink’ term, S(r,t), is present (see Eq. (3)), the 
solution is in the form 


p(r,t) = p(r,t) ■ E(r,t) 
where E(r,t) satisfies 
dlnE d\nE 


+ V r 


dt dr 

The solution of (B13) is: 


+ S(r,t) = 0 


(B12) 

(B13) 


In E(r,t) = — / S (H(x + r(r)),t — x)dx. (B14) 

Jo 

Thus, the solution when a sink term is present specifies to: 
p(r, t) = p om (r , t ) ■ e" So (B15) 


This paper has been produced using the Royal Astronomical 
Society/Blackwell Science BT(rX style Hie. 
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